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We calculate the spin transport of hydrogenated graphene using the Landauer-Biittiker formalism 
with a spin-dependent tight-binding Hamiltonian. The advantages of using this method is that it si¬ 
multaneously gives information on sheet resistance and localization length as well as spin relaxation 
length. Furthermore, the Landauer-Biittiker formula can be computed very efficiently using the 
recursive Green’s function technique. Previous theoretical results on spin relaxation time in hydro¬ 
genated graphene have not been in agreement with experiments. Here, we study magnetic defects 
in graphene with randomly aligned magnetic moments, where interference between spin-channels 
is explicitly included. We show that the spin relaxation length and sheet resistance scale nearly 
linearly with the impurity concentration. Moreover, the spin relaxation mechanism in hydrogenated 
graphene is Markovian only near the charge neutrality point or in the highly dilute impurity limit. 

PACS numbers: 72.25.-b, 72.10.Fk, 72.80.Vp 


INTRODUCTION 

Spin transport in graphene has attracted a lot of at¬ 
tention in recent years due to very long spin relaxation 
times and spin relaxation lengths predicted for this ma¬ 
terial mm- The spin relaxation length in graphene has 
been predicted theoretically to be at least 20 /rm [T], 
whereas experimental values are about an order of mag¬ 
nitude lower, typically around 1-4 yim [3H8]. It has been 
ruled out experimentally that this discrepancy is due to 
hyperfine interaction with naturally occurring iso¬ 
tope in graphene [7]. An attempt to explain the dis¬ 
crepancy between theory and experiment based on mag¬ 
netic impurities in graphene has been given by Kochan et 
al. [3] . Magnetic impurities are very common in graphene 
and may for instance be hydrogen adatoms uni , vacan¬ 
cies uniiii] or embedded metal atoms min in graphene 
pores. Kochan et al. find that 0.36 ppm coverage of 
hydrogen adatoms is sufficient to obtain spin relaxation 
times that are in agreement with experiments. Their 
model is based on the Green’s function of a single hydro¬ 
gen adatom in an infinite graphene sheet and multiplying 
their results with the impurity concentration. In effect, 
they do not include interference effects between scatter- 
ers in their model and their model is thus only valid in 
the highly dilute limit. Experimental measurements of 
graphene in the presence of a strong magnetic field con¬ 
firms that the observed low spin relaxation length is, at 
least in part, due to magnetic impurities in graphene |14j . 
Spin transport in hydrogenated graphene was also con¬ 
sidered by Soriano et al. [mils]. Their method is based 
on a mean-field Hubbard Hamiltonian and the real space 
Kubo-transport formalism. They find that a coverage 
of 15 ppm hydrogen adatoms gives the correct order of 
magnitude of the spin relaxation time |16j , which is more 
than an order of magnitude larger than the prediction by 
Kochan et al. Additionally, the functional form of the two 


theoretical results does not resemble the experimental re¬ 
sult. A recent ab initio study of the spin scattering of 
hydrogen adatoms on narrow armchair graphene nanorib¬ 
bons by Wilhelm et al. m has shown that spin scattering 
off a single hydrogen adatom with defect spin oriented 
perpendicular to the electron spin is sufficient to obtain 
spin-flip conductance on the same order of magnitude as 
the spin-conserved conductance. They also showed spin- 
orbit interactions to be negligible compared to exchange 
interactions in the context of spin scattering on hydrogen 
adatoms. 

The spin relaxation length is determined by the decay 
rate of spin polarization. Zurek et al. |18j have found 
through a phenomenological spin interaction Hamilto¬ 
nian that the spin relaxation decay rate depends on the 
distribution of coupling strengths between a spin sys¬ 
tem and an environment with many independent spins. 
In particular, they find that a Gaussian distribution of 
couplings leads to Gaussian decay of the spin polariza¬ 
tion with respect to time, whereas a Lorentzian distri¬ 
bution leads to exponential decay. It is straightforward 
to demonstrate that the spin relaxation of scatterers on 
a classical Markovian chain is also exponential. There¬ 
fore, exponential decay of spin polarization is typically 
referred to as Markovian behavior [TH]. 

In this paper, we calculate the spin-dependent elec¬ 
tron transport on graphene with hydrogen adatoms using 
the Landauer-Biittiker formalism, which is a widely used 
method for calculating quantum transport in nanoscale 
devices [IIlEoHSl!. We use hydrogen adatoms as they 
are very common magnetic defects on graphene. They 
have a local magnetic moment of approximately 1 ^lb per 
adatom. Additionally, due to local sp^ hybridization, hy¬ 
drogenated graphene has an energy gap |5S]. In particu¬ 
lar, we will demonstrate that the Landauer-Biittiker for¬ 
malism can be used to extract the spin relaxation length 
of a system. We will demonstrate that the spin relaxation 
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is not always Markovian and that spin relaxation length 
and sheet resistance scales nearly linearly with impurity 
concentration. 

THEORETICAL METHODS 

We employ a third-nearest neighbor 7r-electron tight- 
binding (TB) model to set up a spin-dependent system 
Hamiltonian on the form 
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i.o- 

where {cj,<j) is the creation (annihilation) operator 
on the lattice site i (j) with spin cr, is the spin- 
dependent hopping parameter between site i and j and 
Si is the spin-dependent on-site energy on site i. The 
notation (i,j) is used to represent up to third nearest 
neighbor indices. TB parameters of hydrogen adatoms on 
graphene are obtained by fitting the TB band structure 
and local density of states (LDOS) to match the DFT 
band structure and atom-projected partial DOS (pDOS), 
respectively. 

The DFT calculations are carried out using the FHI- 
aims package [55]. It is an all-electron code with nu¬ 
merical atom-centered basis functions. We use the de¬ 
fault tight basis set for each atom type in a spin- 
polarized calculation. The electron-electron interactions 
are treated at the level of the Perdew-Burke-Ernzerhof 
(PBE) exchange-correlation functional |57|. The hydro¬ 
gen adatom on graphene is relaxed in a supercell with 
2x8x8 = 128 carbon atoms, until the forces between 
atoms are smaller than 10“^ eV/A. We expect this su¬ 
percell to be large enough for finite size effects to be 
negligible. Moreover, the DFT self-consistency cycle is 
considered converged if, among other things, the total 
energy changes by less than 10“® eV. We use an 8 x 8 x 1 
A:-point Monkhorst-Pack grid under relaxation. The fi¬ 
nal band structure and pDOS calculations are computed 
using k-grids of 15 x 15 x 1 and 12 x 12 x 1 /c-points, re¬ 
spectively. For the band structure fit, we compare the six 
lowest unoccupied and six highest occupied bands and fit 
the two-dimensional band energies in the first Brillouin 
zone. 

The defect spin moment may not be aligned with 
the chosen quantization axis of the electron spin in 
the leads. In order to take this into account in our 
model, we rotate the defect spins in spin space. To 
do so, we start by writing the spin-dependent Hamil¬ 
tonian as a sum of carbon He and defect Hdefect parts, 

H = Egraphene ^c+Edefects ^defect' The defect Hamil¬ 
tonian is separated into spin channels and is written as 

^defect t/t = Hi ± Ai. By rotating the defect spins in¬ 
dividually on a Bloch sphere with polar angles Oi and 


FIG. 1. Spin-dependent transport is equivalent to having two 
separate channels that couple only at magnetic impurity sites. 


azimuthal angles (f>i it is straightforward to demonstrate 
that the defect Hamiltonians become 


Hdelects ~ Hi + AiiSl Ri, 


cos 9i e sin 9i 

gfoi gjj^ 0. _ pog 0’ 


( 2 ) 

where Hi = (^defect,t+ ^defect.;)/2 is the mean Hamilto¬ 
nian of the spin channels and 2Ai = “ ^defect,; 

is the difference. This way of rotating the defect spins is 
equivalent to the method used in Ref. except that 
this explicitly allows us to rotate each defect spin in¬ 
dependently of the others. It follows from Eq. that 
spin flipping only occurs when the electron spin is not 
aligned with the defect spin, and the coupling between 
spin-channels is proportional to A^. The spin-dependent 
transport is thus equivalent to having two separate chan¬ 
nels that couple when A is non-zero, such as at magnetic 
impurity sites as illustrated in Fig. 

The transmittance between any two leads p and q 
of a multi-terminal system can be calculated using the 
Landauer-Biittiker formula [55] Tpq = TrlFpGFqG'^}, 
where G = [{E + ie)I — H — J2n is the retarded 

Green’s function. and r„ are the self-energy and 
linewidth functions, respectively, of lead n. A small imag¬ 
inary part e = 10“® eV is added to the energy for nu¬ 
merical stability. If the spins are decoupled in the leads, 
it is easy to demonstrate that the spin-channel resolved 
transmittance between the leads of a spin-dependent two- 
terminal system becomes [niEH] 


=Tr{TH^Gri^'>G^}, (3) 

where F^^^ is the linewidth function of the left 

(right) lead with spin a (cr'). The transmittance on 
this form can be computed efficiently using the recur¬ 
sive Green’s functions (RGF) technique (as outlined in 
Refs. [53[ [29]). All calculations are performed on unit 
cells with a relatively large width of 12.8 nm in order 
to minimize finite size effects. Furthermore, the calcula¬ 
tions are performed using periodic boundary conditions 
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transverse to the transport direction and the results are 
averaged over 29 fc-points. Exactly at the CNP, the only 
propagating mode in the leads is at fc = 0. It is therefore 
important to ensure that this is included. 

The spin-conserved transport is defined as T^c = + 

and the spin-flipped transport is defined as Tgf = 
We expect the total transport T = Tgc + Tsf 
to be either Ohmic or localized. For Ohmic transport the 
resistance is R{L) = Rc + R^L/W, where Rc is the con¬ 
tact resistance, Rg is the sheet resistance, L is the device 
length and W is the width of the unit cell. In the lo¬ 
calization regime, the resistance is R{L) = i?cexp(L/^), 
where ^ is the localization length. By fitting the total 
transport to a compound expression 

R{L) = ^ = i?e exp(L/e) + RsL/W , (4) 

we obtain both localization length and Ohmic resistance. 
In the limits ^ oo and Rg —>■ 0, this expression reduces 
to the Ohmic and localization regimes, respectively. 

We can use the spin polarization P to obtain the spin 
relaxation length Ag. According to Zurek et al. [18], 
the spin relaxation mechanism can be either exponen¬ 
tial or Gaussian, depending on the distribution of spin- 
couplings to an environment. In order to include both 
cases as well as any intermediate relaxation mechanism, 
we fit the spin polarization according to the following 
expression 


p(l) ^ '^scjL) Tgf{L) ^ 
^ ^ Tg,{L)+Tgf{L) 


(5) 


It follows that the spin relaxation behavior is exponential 
when n = 1 and Gaussian when n = 2. 


RESULTS 

By fitting to the pristine graphene DFT band struc¬ 
ture, we find the G-G hopping parameters ti = —2.855 
eV, t 2 = —0.185 eV and ts = —0.190 eV for the first, 
second and third nearest neighbors, respectively. The ti 
and t 2 parameters are fitted freely, and ts is included by 
assuming that = ti(0.18/2.7), where the factor is mo¬ 
tivated by earlier models m- The C on-site energy is 
vanishing. By fitting to the band structure of a system 
with an H adatom, the nearest-neighbor C-H hopping pa¬ 
rameter can be taken as spin-independent with a value 
of 9.475 eV. The on-site energy at the hydrogen site is 
spin-dependent with a value of 4.689 eV for the major¬ 
ity spin component and 1.853 eV for the minority spin 
component. As the only spin-dependent parameter is on 
H adatoms, spin-flipping only occur at these sites in our 
model according to Eq. [2|The fitted band structure and 
LDOS are shown in Fig/M where they are compared to 


the DFT calculations. The figure shows excellent agree¬ 
ment between TB and DFT band structures as well as 
the TB LDOS and the DFT pDOS. 



FIG. 2. Spin-dependent band structnre (left) and local den¬ 
sity of states on the H atom (right) for an 8 x 8 graphene 
nnit cell with one H adatom. Solid lines are DFT results and 
dashed lines are TB results. 


The spin polarization of a system containing a single 
H adatom is shown in Fig. When there is only a sin¬ 
gle defect, the transport properties do not depend on the 
azimuthal defect spin angle (j). Therefore, only the polar 
angle 9 and the energy E are varied. The figure shows 
that the spin scatters very strongly near the charge neu¬ 
trality point (GNP), £1 ~ 0, resulting in a significantly 
decreased spin polarization. This is a consequence of 
scattering on H adatoms, which as defect bands that span 
approximately ±0.1 eV around the GNP, cf. Fig. This 
means that a single H adatom with in-plane defect spin 
is able to destroy almost half of the spin polarization for 
energies near the GNP. This is in good agreement with 
Wilhelm et al. im, who found that an = 11 arm¬ 
chair graphene nanoribbon with a single H adatom with 
in-plane spin can have spin-flip transmittance that can 
surpass the spin-conserved part. 

In order to obtain information on the interference ef¬ 
fects of spin-flipping, we calculate the spin polarization 
of a system with two H adatoms separated by a distance 
of 2.21 nm parallel to the transport direction, see Fig.|^ 
The spin polarization is evaluated at the GNP and the 
orientation of the defect spins have been chosen to be 
(01, (/i) = (90°, 0°) and (02,'/ 2 ), respectively. The figure 
shows that the total spin polarization is minimal when 
the defect spins are in-plane and point in the same di¬ 
rection, whereas it is maximal, when the spins are in¬ 
plane and point in opposite directions. In Fig. we saw 
that a single defect with in-plane spin could flip almost 
half of the electron spin to the opposite channel. Now 
we see that by having two defects with oppositely ori¬ 
ented in-plane spins, the second can almost completely 
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FIG. 3. Spin polarization as a function of energy and angle of 
a graphene system with a single H adatom. The inset shows 
an illustration of the device. 
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FIG. 4. Spin polarization of a system with two H adatoms 
with defect spin angles (0i,^i) = (90°,0°) and (02,^2), re¬ 
spectively, at an energy of E = 0.0 eV. The H adatoms are 
placed on a line parallel to the transport direction 2.21 nm 
apart. The inset shows an illustration of the device. 


negate the first spin flip. When the two defect spins point 
in opposite directions, the phase change associated with 
spin-flipping will have equal size and opposite sign. This 
means that the electron spin will be in phase with the 
input spin after the second spin-flip, leading to construc¬ 
tive interference. This is not the case when the defect 
spins point in the same direction. The interference be¬ 
tween defects is thus very important and should not be 
ignored. 

We now turn to calculating the effects of multiple mag¬ 
netic hydrogen adatoms on graphene. We place hydro¬ 
gen adatoms at random positions uniformly distributed 
across the device according to a predefined impurity con¬ 
centration rj. The impurity concentration is related to 
the impurity density by n = 4r7/(-\/3a^), where a is the 
graphene lattice constant. We wish to keep the device 


non-magnetic in order to isolate spin relaxation from 
other magnetic effects. Therefore, we choose the direc¬ 
tions of the defect spins at random, uniformly distributed 
on a Bloch sphere. The transport is calculated for very 
long devices of 147.5 nm, which contain a total of 72,000 
carbon atoms in the unit cell. Using the RGF method, 
we can extract the transport after each slice of the device, 
allowing us to obtain the length-dependent transport di¬ 
rectly. In order to minimize the effects of the finite width 
of the unit cell, we average over an ensemble of 150 de¬ 
vice realizations. The spin polarization as a function of 
length and energy for different impurity concentrations 
is shown in Fig.j^as well as an example of transmittance 
and spin polarization as a function of length for a single 
energy and impurity concentration. We show the loga¬ 
rithm of the spin polarization in the range between -1 and 
0 in order to highlight the spin relaxation length, which 
is defined as the length at which ln(P(L)) = —1. The fig¬ 
ure shows that the spin polarization decays very fast for 
energies close to the H adatom defect bands, cf. Fig. 

As expected, the spin polarization decays faster with in¬ 
creasing impurity concentration. Note that the spin po¬ 
larization also decays for energies away from the H defect 
bands, due to the relatively small spin splitting in the 
remaining band structure. The small energy-dependent 
oscillations in Figs. and are due to finite size effects 
originating from the finite width of the unit cell. 

By fitting the spin polarization to Eq. [^we obtain the 
spin relaxation length as well as the exponent n, which 
provides information on the spin relaxation mechanism, 
see Fig. Hi- A few examples of the fitting procedure is in¬ 
cluded in Fig. in order to illustrate the excellent qual¬ 
ity of the fits. The carrier concentration on the figure is 
computed at Fermi energies corresponding to the energy 
axis. Positive carrier densities refer to electron doping 
and negative carrier densities refer to hole doping. The 
spin relaxation length is very short for energies near the H 
defect bands. For the same energies, the spin relaxation 
mechanism is predominantly exponential with an expo¬ 
nent of n ~ 1. For energies further away from the CNP, 
the spin relaxation length increases. We note that As has 
two minima near the CNP, which are correlated with the 
large spin-splitting of the H adatom defect bands. Ex¬ 
actly at the CNP, the spin-splitting of the defect bands is 
vanishing, resulting in a local maxima. The figure shows 
that there is an almost linear scaling of the spin relax¬ 
ation length with respect to impurity concentration, es¬ 
pecially near the CNP. Away from the CNP we observe 
that n decreases with decreasing impurity concentration. 
This suggests that the spin relaxation mechanism tends 
toward exponential (Markovian) behavior in the highly 
dilute impurity limit. Importantly, we see that the decay 
of the spin polarization as a function of length need not 
be exponential nor Caussian, which means that any good 
theory on spin relaxation should not presume anything 
about the spin relaxation behavior, except in the limit 
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FIG. 5. Ensemble-averaged transmittance and spin polariza¬ 
tion as a fnnction of length (top) for a system with impurity 
concentration t] = 500 ppm calculated at the CNP. Ensemble- 
averaged spin polarization as a function of energy and length 
(bottom) for different impurity concentrations. The dashed 
lines show the spin relax;ation length. 


of very dilute systems, where the approximation of ex¬ 
ponential decay seems to be valid. For energies near the 
CNP, the normalized localization length is Asry « 0.01 
nm. In order to obtain experimentally observed spin 
relaxation length of about Xs — 2 /rm [3], the impu¬ 
rity concentration should be ry « 5 ppm, which is more 
than an order of magnitude larger than the prediction by 
Kochan et al. [3] of 0.36 ppm. We expect our model to be 
more accurate as it is based on a full transport calcula¬ 
tion and therefore takes interference effects into account. 
Our prediction of the impurity concentration is, however, 
in closer agreement with Soriano et al. |16j . who found 
that an impurity concentration of 15 ppm gives spin re¬ 
laxation times in agreement with experiment based on 
time propagation of the spin polarization operator using 
a self-consistent Hubbard model. 

A comparison of spin relaxation lengths obtained by 
the current model and those obtained by other theoret¬ 
ical methods and experiment is presented in Fig. We 
have rescaled our 500 ppm result to 5 ppm by multiplying 
it by a factor of 100. The two other theoretical results 
mm have calculated the spin relaxation time ts, which 
is related to the spin relaxation length by Xs = vsts in 
the ballistic regime and by Ag = ^/Dsts in the diffu¬ 
sive regime, where vs is spin carrier velocity and and Ds 


FIG. 6. (a) Normalized spin relaxation length Xsr^ and expo¬ 
nent n (inset) obtained by fitting against Eq. The spin re¬ 
laxation lengths are normalized with the defect concentration 
in order to illustrate their nearly linear scaling with respect 
to it. (b) Examples of fitting the spin relaxation against Eq. 
of a system with 5000 ppm H adatoms for different energies. 
The energies are between 0.0 eV (fastest decay) and 1.0 eV 
(slowest decay) in steps of 0.2 eV. The dots are the ensemble 
averaged spin polarizations and the lines are the correspond¬ 
ing fitted functions. 


is the spin diffusion constant. In the low-defect-density 
case, we expect to be in the ballistic regime. Therefore we 
compare results that are all obtained in the low-defect- 
density case. We obtain a velocity vs = 1.65 x 10^ m/s 
by least-of-squares fitting between our result and the an¬ 
alytic result obtained by Kochan et al. We observe that 
the result by Kochan et al. is in fairly good agreement 
ours regarding the location of the two minima near the 
CNP and in quantitative agreement further away from 
the CNP. However, their result predicts variations over 
several orders of magnitude near the CNP, whereas our 
result predicts a variation of only about a factor of 2. In 
fact, their result is singular exactly at the CNP. Further¬ 
more, we compare with experimental results of hydro¬ 
genated graphene obtained by Wojtaszek et al. [5]. Note 
that the experimental results were obtained without de¬ 
tailed knowledge of the defect concentration. However, 
the authors estimated the concentration to be 200 ppm. 
Lastly, we compare our results to the theoretical result 
by Soriano et al. [HI- The figure shows that their re¬ 
sult is neither in qualitative agreement with our model 
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FIG. 7. Comparison of spin relaxation lengths obtained by 
different authors. 

nor the analytic result by Kochan et al. or experiment. 
We speculate that the deviation arises from the fact that 
Soriano et al. uses vacancies in graphene to model hy¬ 
drogen adatoms, whereas both our model and the model 
used by Kochan et al. employ a parametrization of hy¬ 
drogen on graphene. Theoretical predictions [SIlelEI] 
including our own, show that the spin relaxation time 
(or spin relaxation length) decreases with increasing im¬ 
purity concentration. However, experimental work on 
hydrogenated graphene shows that the spin relaxation 
time (or spin relaxation length) actually increases with 
increasing impurity concentration [3]. The origin of this 
discrepancy remains elusive, but could stem from inter¬ 
action between graphene and the substrate, as this has 
not been included in any of the theoretical models. 

By fitting the total transmittance against Eq. we 
obtain the Ohmic sheet resistance as well as the local¬ 
ization length, see Fig. We observe localization near 
the H defect bands, cf. Fig. and vanishing localiza¬ 
tion elsewhere. Additionally, the figure shows that the 
sheet resistance scales linearly with respect to impurity 
concentration. However, the scaling of the localization 
length is far from linear, which shows that the localiza¬ 
tion induced per atom decreases with increasing impurity 
concentration. Furthermore, as the impurity concentra¬ 
tion is decreased the energy window, at which there is 
localization narrows. 


CONCLUSIONS 

We have demonstrated that the Landauer-Biittiker for¬ 
malism can be used to calculate spin-dependent transport 
of systems with magnetic impurities with individually 
oriented magnetic moments. In this work, we study hy¬ 


Carrier density [xlO^^/cm^] 



Energy [eV] 


FIG. 8. Normalized Ohmic sheet resistance Rs/rj (top) and 
normalized localization length ^/rj (bottom) obtained by fit¬ 
ting against Eq. ^ 


drogen adatoms on graphene. By calculating the length- 
dependent transport, we can extract properties such as 
spin relaxation length, localization length and sheet re¬ 
sistance. We have shown that there is strong local¬ 
ization for energies around the hydrogen-induced defect 
bands, which also leads to a very high sheet resistance. 
Away from the defect bands there is vanishing localiza¬ 
tion. Furthermore, we have demonstrated that the spin 
relaxation length is very short for energies around the 
hydrogen-induced defect bands and that the spin relax¬ 
ation mechanism is exponential (Markovian) near the 
CNP and non-exponential (non-Markovian) otherwise. 
Additionally, we have shown that spin relaxation length 
and sheet resistance scale nearly linearly with impurity 
concentration, whereas the localization length does not. 
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